Non-Markovian quantum state diffusion for absorption spectra of molecular 

aggregates 



O 

(N 

< 

o 

(N 



a,: 
c ■ 

ctf ■ 
=5 

cr 



> 

\o 
in 
cn 

od 

o 
o 



X 



Jan Rodcn, 1 Walter T. Strunz, 2 and Alexander Eisfeld^B 

1 Max- Planck- Institut fur Physik kornplexer Systeme, Nothnitzer Str. 38, D-01187 Dresden, Germany 
Institut fur Theoretische Physik, Techmsche Universitat Dresden, D-01062 Dresden, Germany 

(Dated: August 24, 2010) 

In many molecular systems one encounters the situation where electronic excitations couple to 
a quasi-continuum of phonon modes. That continuum may be highly structured e.g. due to some 
weakly damped high frequency modes. To handle such a situation, an approach combining the non- 
Markovian quantum state diffusion (NMQSD) description of open quantum systems with an efficient 
but abstract approximation was recently applied to calculate energy transfer and absorption spectra 
of molecular aggregates [Roden, Eisfeld, Wolff, Strunz, PRL 103 (2009) 058301]. To explore the 
validity of the used approximation for such complicated systems, in the present work we compare 
the calculated (approximative) absorption spectra with exact results. These are obtained from the 
method of pseudomodes, which we show to be capable of determining the exact spectra for small 
aggregates and a few pseudomodes. It turns out that in the cases considered, the results of the two 
approaches mostly agree quite well. The advantages and disadvantages of the two approaches are 
discussed. 



I. INTRODUCTION 

There is growing interest in systems composed of indi- 
vidual monomers that interact via resonant dipole-dipolc 
interaction. Upon electronic excitation this transition 
dipole-dipolc interaction between the monomers is re- 
sponsible for a collective behaviour of these systems. 
Besides the classical examples like Van-der-Waals crys- 
tals [l|, , aggregates of organic dyes , and light har- 
vesting units of plants, algae, and bacteria @,|I|-[Tl[ many 
new systems have emerged. Examples are ultra-cold Ry- 
dberg atoms [l2l - [l6l quantum dots [l?], HH , assemblies 
of na no-p articles [IS . |20| . and recently also hybrid sys- 
tems mf. 

The common approach to describe these aggregates is 
to treat the monomers as electronic two-level systems. 
Besides the electronic degrees of freedom, however, one 
often has to take into account nuclear degrees of freedom 
explicitly. For instance in the case of molecular aggre- 
gates 0, |22M2(| , which will serve as the primary exam- 
ple in this work, the electronic excitation of a monomer 
couples strongly to internal vibrational modes of the 
monomer and to modes of the surroundings. 

Often, the monomer spectrum is dominated by one vi- 
brational progression which is considerably broadened. 
A commonly applied approximation is then to only con- 
sider one effective mode corresponding to this progres- 
sion [U HH and to take the broadening into account by 
convoluting with some lineshapc function which is usu- 
ally assumed to be Gaussian. It has been shown that 
using this approach already important features of exper- 
imental spectra can be reproduced [28H331 ]. Although 
the resulting spectra reveal many characteristics of the 
aggregates, important aspects like the detailed shapes of 
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the J-band [34[ and the H-band [35[ cannot adequately 
be described by considering only one vibrational mode. 
On the other hand, the exact inclusion of only one vibra- 
tional mode already complicates the treatment of molec- 
ular aggregates considerably, so that this approach is re- 
stricted to small aggregates. This problem becomes even 
more serious when one attempts to include more modes 
in this manner. 

When the interaction of the chromophorc monomers 
with the environment is negligible (e.g. in high resolu- 
tion spectroscopy in helium nanodroplets [24l. |36|). then 
the explicit inclusion of vibrational modes is of great im- 
portance. However, for typical spectra in solution or in a 
solid state matrix (where a strong coupling between the 
chromophores and the environment is present) it seems 
better to use a continuum of vibrations that couple to 
the electronic excitation to account for the large number 
of environmental degrees of freedom. 

This interaction between the electronic excitation and 
the vibrations is conveniently encoded in the so-called 
spectral density. It describes the frequency-dependent 
coupling between the system (the electronic degrees of 
freedom) and the (continuum of) harmonic oscillators. 
In the Markov case the spectral density is assumed to 
be flat in the relevant frequency regions. Clearly, for 
the considered monomers this assumption does not hold. 
Due to strong interaction with some internal vibrational 
modes, the spectral density will be highly structured (i.e. 
frequency-dependent), indicating that a non-Markovian 
theoretical framework is required. 

An approach to tackle this complicated problem was 
recently presented in Ref [37|. The method is based on 
the non-Markovian quantum state diffusion (NMQSD) 
description of open quantum systems [38| ■ Here, the sys- 
tem part is chosen to contain only the electronic degrees 
of freedom which interact with a non-Markovian envi- 
ronment (the bath) comprising all vibrations (internal 
modes of the monomers as well as external modes). One 
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then can derive a stochastic evolution equation for states 
in the (small) space of the system part. However, solv- 
ing the exact evolution equation turns out to be very 
difficult due to the appearance of a functional derivative 
w.r.t. functionals containing the bath degrees of freedom. 
To overcome these difficulties, in an approximation only 
the zeroth order of a functional expansion (we will refer 
to it as ZOFE approximation) of the problematic term 
is taken into account [13, [39[ . For several (simple) prob- 
lems this procedure has been shown to give the exact 
result [13,|4l|. However, for more complex problems like 
the molecular aggregates studied in this work, the range 
of validity of the approximation is not clear. It should be 
noted that the NMQSD approach in combination with 
the ZOFE approximation provides a very efficient calcu- 
lation scheme: in order to obtain the absorption spec- 
trum of the aggregate and the energy transfer between 
the monomers, the equations one has to solve are in the 
small Hilbert space of the electronic degrees of freedom 
alone [37j . 

One aim of the present paper is to examine the valid- 
ity of the ZOFE approximation leading to the calculation 
scheme presented in Ref. [37[. To this end, we compare 
with an approach where so-called pseudomodes |42h44| 
are included into the system part together with the elec- 
tronic degrees of freedom. The electronic degrees of free- 
dom now couple only to the pseudomodes, the pseudo- 
modes in turn are then coupled to a Markovian bath. For 
a spectral density consisting of a sum of Lorentzians the 
pscudomode method is exact (taking one pseudomode 
for each Lorentzian). This allows to directly compare 
the approximative NMQSD-ZOFE treatment with exact 
calculations. However, due to the inclusion of the pseu- 
domodes into the system part, the numerical solution 
of the corresponding evolution equation is limited to a 
rather small number of monomers in the aggregate with 
only a few pseudomodes, i.e. only a few Lorentzians in 
the spectral density. 

Besides the possibility of comparing the NMQSD- 
ZOFE approach with exact calculations, the pseudomode 
method has also some physical significance: one can think 
of the pseudomodes as internal vibrational modes of a 
chromophore that strongly couple to the electronic ex- 
citation and which are damped by the coupling to the 
remaining vibrations. 

The comparison between zero temperature absorption 
spectra of small aggregates calculated using the NMQSD- 
ZOFE approach and spectra obtained from the exact 
pscudomode approach shows that in the cases consid- 
ered there is mostly quite good agreement between the 
two approaches. We will discuss in which situations the 
approximative result of the NMQSD-ZOFE approach is 
expected to deviate from the exact solution. 

The structure of this paper is as follows: In Section [TTl 
we introduce the Hamiltonian of the aggregate. The 
Hamiltonian is written as the sum of a system part (con- 
taining only electronic degrees of freedom), an environ- 
mental part (containing all vibrational modes), and the 



part of the interaction between electronic degrees of free- 
dom and vibrations. In the following Section IIII1 the 
basic formulas that are used to calculate the absorp- 
tion spectrum are given by specifying the initial state 
and introducing the dipole correlation function. In Sec- 
tion IIV1 the general Non-Markovian Schrodingcr Equa- 
tion (NMQSD) approach is applied to the case of an 
aggregate. It is shown how the absorption spectrum 
can be obtained in this approach. Next, in Section [V] 
the ZOFE approximation is introduced. Then, in Sec- 
tion IVI1 the pseudomode (PM) approach is presented. 
In Section IVTI1 the NMQSD-ZOFE absorption spectra 
are compared with exact PM spectra. We conclude in 
Section [Villi by summarizing our findings. Details of the 
calculations and minor results have been placed in the 
appendices. In Appendix \^ two exactly solvable cases 
(namely that of non-interacting monomers and the case 
where the coupling to the vibrations can be considered to 
be Markovian) are discussed. In Appendix [B] it is shown 
how to obtain the absorption spectrum using the PM 
approach. The numerical implementation is discussed. 
Finally, in Appendix [Cl it is shown that for a Lorentzian 
spectral density the absorption spectrum obtained from 
the exact NMQSD approach is equal to the spectrum 
obtained from the PM method. 



II. THE AGGREGATE HAMILTONIAN 

We consider an aggregate consisting of N monomers, 
labelled by n = 1, . . . , N. For each monomer n we take 
into account its electronic ground state | <\fl n } and one 
excited electronic state | 4>n ) • The transition energy be- 
tween these two states (whose wave functions we take 
to be real) is denoted by e n . Apart from the two elec- 
tronic states, each monomer has a collection of vibra- 
tional modes comprising internal modes as well as modes 
of the local environment of the monomer. We will re- 
fer to these degrees of freedom as "nuclear coordinates" . 
The electronic excitation of monomer n couples to its 
vibrations and we choose the vibrational modes to be 
harmonic and the couplin g to be linear (making contact 
to previous work [H, IMTMSl ) • The Hamiltonian of 
monomer n is then given by 



(i) 



with the Hamiltonian of the vibrations in the electronic 
ground state 



H 9 n = ^2 twn\al x a n \ 



(2) 



(the energies of the vibrational ground states of all modes 
in the electronic ground state are chosen to be zero) . The 
Hamiltonian of the vibrations in the excited electronic 
state comprising a shift, reads 



Hi 



A 



TuQn\a\ lX a n \ - ^2 K «A(ai A + a n \). (3) 
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Here a n \ denotes the annihilation operator of mode A of 
monomer n with frequency cu n \. The coupling strength 
with which electronic excitation of monomer n couples 
to mode A of this monomer is denoted by K n \- 

For the aggregate we assume that the electronic wave 
functions of the monomers do not overlap. The electronic 
ground state of the aggregate is then taken as the product 



N 



isei>= n 



(4) 



m— 1 



of the electronic ground states | 4> m } of all monomers. 
A state of the aggregate in which only monomer n is 
electronically excited and all other monomers are in their 
electronic ground state we denote by 



N 



(5) 



We expand the aggregate Hamiltonian w.r.t. the states 
Eqs. (j4]) and (J5j and neglect states with more than one 
electronic excitation on the aggregate. Thus we obtain 
the Hamiltonian 



JV 



H = H 9 \g cl )(g cl \+H e J2 1 7r„)(7r„ 



(6) 



for the aggregate, with the part 

N 



(7) 



and a part H env describing the "environment" of vibra- 
tional modes 



(11) 



ri=l A 



The coupling of electronic excitation to these vibrations 
is expressed through 



N 



H in t = -y^Kn)(7Tniy^ ^nX{a n \ + 0>n\)- (12) 



We emphasize that with Eqs. (|10| -(|T2 j) we make a special 
choice of the three parts, system, interaction and envi- 
ronment, of the aggregate Hamiltonian. In particular the 
system part Eq. (fTOf contains only electronic degrees of 
freedom. In Scction lVll where we will introduce the pseu- 
domode approach, we will include also vibrational modes 
(pseudomodes) into the system part. 

A useful quantity describing many aspects of the cou- 
pling of the system degrees of freedom to the vibrational 
environment is the so-called spectral density [BlT ]. which 
for monomer n is given by 



J n (uj) = ^ \ K nx\ 2 <5( w ~ Unx) 



(13) 



and which will be used later, generalized to a continuum 
of vibrational modes. 



III. ABSORPTION OF THE AGGREGATE 



for the electronic ground state and the part 

N / N \ 

ffe =E \ H n + E H rn) M(*n 

N 

+ E ^ nm \' K n){ Km 
n,m—l 



(8) 



for the electronically excited state. The matrix element 
V nm , causing electronic excitation to be transfered from 
monomer n to monomer m via transition dipolc-dipolc 
interaction, is taken to be independent of nuclear coordi- 
nates (note that H 9 and depend on nuclear coordi- 
nates through Eqs. © and ©). With the Hamiltonians 
of the monomers Eqs. (|2|) and Q we can write Eq. ([§]) 



H e — H sys + H int + H env , 



with the purely electronic "system" part 

N N 



(9) 



H sys = ^ E n \ TT„ )(7T n I + ^ V nm \-K n ){-K m |, (10) 



We consider absorption of light by the aggregate at 
zero temperature. Initially, the aggregate is taken to be 
in its total ground state 



|*(t = 0)) = |ffel>|0), 



(14) 



which is a product of the electronic ground state | <7 e i ) 
defined in Eq. (|4]) and the ground state | ) = Yl nX | 0„a ) 
of -ffenv , i-e. all vibrational modes of all monomers are in 
their ground state | 0„a ) ■ The absorption of light with 



frequency v and polarization £ is then given by [51 
A(v 

Here, 



Re / dt e wt M{t). 



(15) 



M(t) = ($(t = 0)\jl-£ e- lHt / h jt-E*\${t = Q)) (16) 

is the dipole correlation function where /2 denotes the 
total dipole operator of the aggregate, given by the sum 



n=l 



n,m— 1 



N 
n=l 



(17) 
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of all monomer dipole operators jt n . Note that the Hamil- 
tonian H in the exponent of the propagator in Eq. (|16l) is 

the aggregate Hamiltonian of Eq. (|6|) . We take ft n to be 
independent of nuclear coordinates. Then the correlation 
function Eq. fll6|) with Eq. ([14]) reads 

N 

M(t) = OC ■ ?){ Km I ( I e^ m ' h | ) | n n ) (ft, • £*), 

n,m— 1 

(18) 



with the monomer transition dipoles 



fl n = {<jf n \fln\K)- 



Defining 



Mtot 



JV 



y ] iMn ' & 
\ ri=l 



* 1 2 



(19) 



(20) 



equation (|18[) can be written as 

M{t)=n 2 tot (*{t = Q)\^{t)) (21) 

with 

) = cxp(-iH e t/h)\ ^>(t = 0) ), (22) 
where H e is given by Eq. © and the initial state is 

|*(t = 0)> = |^)|0>. (23) 



The normalized initial electronic state | V'o ) in Eq. ([23 
is given by 



JV 



I V>0 ) = — VVn ' £*)\ 7T„ } 
Mtot , 



(24) 



and contains explicitly the geometry of the aggregate 
via the orientation of the transition dipoles fl n of the 
monomers. In the following, our goal will be to obtain 
the state \^(t)) to be able to calculate the absorption 
spectrum according to Eqs. ([T5|) and ([2"Tj) . 



IV. THE GENERAL NMQSD APPROACH 

The correlation function M{t) of Eq. ([2"Tjl can be ob- 
tained using the framework of stochastic Schrodingcr 
equations, here the non-Markovian quantum state dif- 
fusion (NMQSD) approach (Ref. [H). Note, however, 
that no stochasticity will enter in the following calcula- 
tions due to the fact that for a zero temperature absorp- 
tion spectrum considered in this work, one has to project 
on the environmental ground state (see section IIVBI) . 
For energy transfer between monomers, however, the full 
stochasticity of the non-Markovian quantum state diffu- 
sion (NMQSD) approach (Ref. [38|) resurfaces (37I . 



A. NMQSD evolution equation for system states 

We briefly summarize the NMQSD approach as applied 
to molecular aggregates. First, we transform H e to the 
interaction representation w.r.t. H cnv , yielding 

N 

H e (t) - H sys + ]T (i„4 (t) + L\A n (t)) (25) 



n = l 



with 



and 



L n = - | TTn } ( 7T„ I = L n 



^nX^nX^ 



(26) 



(27) 



Thus, we have the time-dependent Schrodingcr equation 

ihdtW)) =H e (t)\V(t)) (28) 

in the interaction picture. In the next step we expand 
the total wave function | ) w.r.t. Bargmann coherent 
states | z n \ } = exp(z n Aa^ A )| n \ ) [52[ of the vibrations, 
where | 0„a ) is the ground state of vibrational mode A of 
monomer n and the z n \ are complex numbers. One then 
obtains 



*(*)> 



d 2 z 



(29) 



with states \ifj(t,z*)) in the space of the electronic 
system H sys , |z) = H n ]J X \ z n \), and d 2 z = 
dRe(z) dlm(z). Inserting Eq. (jiSJ into Eq. (gSJ) one 
finds [38| that the states | ip(t, z*) ) appearing in Eq. ([29]) 
obey an evolution equation 



W(f,s*)) 



H sys \i>(t,z*))+J2L n zlMt,z*)) 



1 v^ rt , , J\tp(t,z*)' 



(30) 



in the small Hilbcrt space of the electronic degrees of 
freedom alone, with time-dependent complex numbers 



iy 



(31) 



and where 

a n (t-s) = (A n (t)At(s)) T=0 = (0\A n (t)Al(s)\0) 

(32) 

is the zero temperature bath correlation function of 
monomer n, which encodes the interaction of an elec- 
tronic excitation with the environment of vibrations. 
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Note that in this case of zero temperature, a n (t—s) is just 
the Fourier transform of the spectral density Eq. (|13p . i.e. 
a„(r) = Jdu e' WT J n {uj). 

For a non-Markovian bath, the solution of Eq. (|3"U| is 
complicated due to the appearance of the memory inte- 
gral over a n (t — s) and the appearance of a functional 
derivative as an integrand. To deal with the functional 
derivative we follow Rcf. [40(1 and write 



Sz 



(33) 



with an operator acting in the space of the electronic 
system part H sys only. The operators OW(i,s,z') have 
to obey the consistency condition 



(oW(t, S ,z*)|#,z«)))=^a|^,z*)), (34) 



which leads, using Eq. (f3"U)) and introducing the abbrevi- 
ation 

6 (B) (i,z«) = jpf da a n (t-s)0^(t,s,z*), (35) 
to an evolution equation [53j 



dtO (m) (t, s , 



--H sys ,0^ m \t,s,z*) 



E \^nZ* t n 
n 

E L 



iioW(^*),oW( M , z *) 



t 5Q(")(t,z*) 



(36) 



The latter has to be solved with initial condition 

0^(s,s,z*) = L n . (37) 

Finally, using Eq. (|3"3"|) and (|3"5"j) . the evolution equation 
pop turns into the linear non-Markovian QSD equation 



d t \i>(t,z*)) 



-H sys \^(t,z*)) 



Y,(L n zl n -LlO^{t,z*))m,z*)). 



(38) 



Due to the functional derivative appearing in Eq. Q36p . 
the operator 0^(t,z*) Eq. ([33]) cannot be evaluated 
in the general case. However, it can be obtained in 
some special cases, e.g. in the Markovian case (see Ap- 
pendix IA 2p and the case of non-interacting monomers 
(see Appendix I A ip . 



B. Absorption in the NMQSD approach 

Using the NMQSD approach, the correlation function 
M(t) of Eq. (pH]) can be calculated as follows. Inserting 



the expansion Eq. (|29|) into Eq. (|2Tj) yields 
d 2 z 



M(t) = m4 



[ih\W,z*))(0\z), (39) 



where \ipo) is given by Eq. (|24[) . Here we used that 
(0 \e~ lH '* nvt ' ,i - which appears through the transforma- 
tion to the interaction representation - is equal to ( | . 
From equation (|39|) we get [52[ 

M{t) = ^ ot (^{t,z* = 0)> 



where | ip(t, z* = 0) ) can be obtained using Eq. 
z* = 0, i.e. 

d t \^(t,z* = 0)) = 



(40) 
with 



y- n Hsys - E i » (5(n) (*> z * = °)J i ^ z * = °) ) 

(41) 

with the initial condition | ijj(t = 0, z* = 0) ) = | i/jq )• 
Note that due to the appearance of the functional deriva- 
tive in Eq. (f3"6")l the determination of (5(™)(i) is still a 
formidable task. 



V. THE ZEROTH ORDER FUNCTIONAL 
EXPANSION (ZOFE) APPROXIMATION 

The main task in the NMQSD approach is to obtain 
the operator (5( n )(i,z*) (and 0^(t,s,z*) respectively) 
being defined in Eq. (|35|) (and Eq. (|33|) ). An evolution 
equation for 0^ n '(t, s, z*) is given by Eq. (f36j) . Although 
for the calculation of the zero temperature absorption 
spectrum from Eq. (|4Tj) only the values of 0^ n \t, s, z*) 
for z* = are needed, Eq. (|36p contains the z* depen- 
dence via the functional derivative in a non-local way. To 
simplify Eq. flM]) we follow Ref. [39[ expanding the oper- 
ator 0^(t,s,z*) w.r.t. Zj n ina functional way and keep 
only the zeroth order term of the functional expansion. 
In other words, we approximate 



(Tl) (t,.s,z*) ^0^>(t,s) 



(42) 



to be independent of z* and refer to Eq. (|42j) as the ze- 
roth order functional expansion (ZOFE) approximation. 
Then, from Eq. ( 1361) one obtains an approximate evolu- 
tion equation [39[ 



d t O^(t,s) 



:H sys ,0^(t, S ) 



J2[LlO^(t),0^( t> s) 



with initial condition Oq"^(s, s) = L n (where we obtain 
0^ n \t) from O { n) (t,s) via Eq. (J35J). Inserting the ap- 
proximate operator (t) into Eq. (|4ip gives 
d t \^(t,z* =0)) = 

(-^-E^o n) w) wm* = o)> (44) 



(43) 
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whose numerical implementation is straightforward [67| . 
To obtain the absorption spectrum of the aggre- 
gate, we solve equation (|44|) for the initial condition 
■0(0, z* = 0) ) = | Vo ) given by Eq. flU}. Then, we cal- 
culate the spectrum via Eqs. (|40| and (fT5|) . 



VI. PSEUDOMODE METHOD 

As we will show in subsequent sections, the NMQSD 
approach in combination with the ZOFE approximation 
presented in Section IIVI offers a highly efficient method 
to calculate absorption spectra of molecular aggregates. 
However, it relies on the rather abstract ZOFE approx- 
imation made in Eq. ()42[) . To obtain information about 
the accuracy of the approximation, in the following we 
will compare with exact calculations. It is clear that in 
general, an exact determination of aggregate spectra is 
a formidable task, impossible for arbitrary spectral den- 
sities and an arbitrary number N of monomers. In this 
section we review a method f43 - l44j that allows a numeri- 
cally exact calculation of spectra for small aggregates for 
bath correlation functions of the form 



iQ„j (t— s)— *tnj\t— s\ 



(45) 



which corresponds to the spectral density 



7T ^ J 



7raj 



r '-r' ~ (w — Slnj) 2 + 7 



(46) 



being a sum of Lorentzians centered at fl n j with width 
j n j. For the numerical implementation of the method, 
the number of Lorentzians which can be taken into ac- 
count is limited (see examples in section IVTI)) . As a spe- 
cial case, the limit 7 nj — > is included, i.e. the case of 
undamped vibrational modes, which has been extensively 
studied in the literature [H, S3, H HI . From the 

point of view of open system dynamics, the memory time 
of the environment is clearly infinite in such a case. 

We want to point out that the NMQSD-ZOFE ap- 
proach is not restricted to the special form (|4"5j) for the 
bath correlation function. However, note also that in 
principle an arbitrary bath correlation function can be 
approximated by a sum of exponentials in the form of 
Eq. (gSJ) as discussed in Refs. p3.l43.l57j. 



A. The pseudomode (PM) Hamiltonian 



In the pseudomode approach (see e.g. |42j-|44|]) the sys- 
tem part of the aggregate Hamiltonian is enlarged. Apart 
from the electronic degrees of freedom, auxiliary vibra- 
tional degrees of freedom (pscudomodes) are included in 
the "system" , each coupled to a Markovian bath in a way 
specified below. 



The Hamiltonian of the aggregate is written as 

N 

H = Ha\g cl )(g cl \+H e Y / \Kn)(n n \ (47) 



where H 9 is the vibrational Hamiltonian in the electronic 
ground state. The relevant Hamiltonian H e of the aggre- 
gate in the excited electronic state is given by 



H e = H« 



Ho 



(48) 



with the following choice of system, interaction and en- 
vironment: the system part is chosen to be 



JV 



n=l j 



N 



(49) 



n=l j 



Here, as before, the Hamiltonian H sys is that of the 
purely electronic system given by Eq. (|10p . Additionally, 
for each monomer we include a set of vibrational modes 
(second term of Eq. (|49|0 . where mode j of monomer 
n has a frequency Q, n j (see Eq. (|45|) ) and annihilation 
operator b n j. These modes, enumerated with index j, 
are referred to as pseudomodes (PM). The third term of 
Eq. (|4"!)|) describes the coupling of electronic excitation 
on monomer n to its PM j with coupling strength T n j 
(see Eq. (|4"5|) ). Each of the PMs has its own environment 
whose modes are enumerated with index p, so that in 
obvious notation, 



N 



n=l j p 



njp a njp- 



(50) 



The PMs interact with their environments through 



A' 



E E E {^njpQnjpbnj + 'P®nj p 3 ™! 



(51) 



n=\ j 



where k n j p denotes the coupling strength between PM j 
of monomer n to mode p of its local environment. 

We now take the bath correlation functions of the PMs 
to read 



(*-•) = £ 



2 g -itj, 



=2h 2 lnj 8{t-s), 



(52) 

i.e. the PMs couple to a Markovian environment (the pa- 
rameter 7„j is that of Eq. (J45J ) . For a Markovian envi- 
ronment, however, we are in the regime of the standard 
Markov QSD of Section IIVI and find an exact solution 
(i.e. without applying the ZOFE approximation) for the 
time-dependent aggregate state and the absorption spec- 
trum. Clearly, the price to pay is the need to propagate 
in the much larger Hilbert space of electronic and PM 
degrees of freedom. 
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We show in Appendix [C] that for a bath correlation 
function of the type of Eq. (g5]), the exact NMQSD ap- 
proach without ZOFE approximation of Section IIVI and 
the PM method result in the same absorption spectra. 

B. Absorption spectrum within the PM approach 

As in Eq. (fl"5|) the absorption spectrum is calculated 
from 

/>oo 

A{v) = Re / dt e ivt M(t), (53) 
Jo 

where now instead of Eq. (j40]) we have 

M(t) = tf ot (i> \i,(t)) (54) 

with 

= |V , o)|s , pm)- (55) 

Here, the initial electronic state | ipo ) contains the action 
of the dipole operator, see Eq. ([24| and | <7pm ) is a prod- 
uct of the ground states of all PMs. As shown in detail 
in Appendix [Bl the state | ip(t) } in Eq. (|54|) is obtained 
by solving the Schrodingcr equation 

d t \m > = ( - l -H sya -EE inAjK ) i m > 

V « =1 3 J 

( 56 ) 

with the initial state | ipo ) and -ff sys from Eq. (|4T)|) . We 
describe in Appendix [B] how this differential equation 
is solved numerically. There, we see that its solution 
becomes quite involved for large aggregates due to the 
explicit inclusion of the PMs into the "system" . 

We note that Eq. (|56p is just the Markov QSD equation 
with zero noise, and for the PM method takes the role of 
Eq. gl]) in the NMQSD approach. 

VII. COMPARISON BETWEEN ZOFE AND 
EXACT PM SPECTRA 

In this section, we will compare absorption spectra cal- 
culated using the ZOFE approximation with numerically 
exact calculations obtained in the PM approach for spec- 
tral densities of the form Eq. (|4"oT) . While being exact in 
the PM approach, the price one has to pay here is the in- 
crease of the number of degrees of freedom of the system 
part due to the inclusion of the PM into the system. This 
entails a rapid growth of the Hamiltonian matrix as the 
number 71/ of PMs or the number N of monomers of the 
aggregate is increased. Thus, using the PM approach 
computer capabilities limit us to absorption spectra of 
aggregates with roughly N = 2 and M « 6 or N = 3 
and M s=s 5 etc. The values of N and M we can handle, 
depend also on the coupling strength T n j of the PM to 



the electronic excitation: the larger the coupling T n j , the 
more basis states have to be taken into account. 

In the following, we consider a linear arrangement 
of monomers with identical properties, i.e. T n j = Lj, 
Q n j = ilj, jnj = 7j for the jth PM of all monomers. For 
simplicity, we take all transition dipole moments of the 
monomers to be equal and consider only nearest neigh- 
bour interaction between the monomers, which we de- 
note by V. We follow Simpson and Peterson [58| and 
speak of strong/intermediate/weak interaction, if V is 
larger/similar/smaller than the width of the monomer 
spectrum. As a parameter for the strength of the cou- 
pling of electronic excitation to the PMs we use the di- 
mensionless Huang-Rhys factor [59[ 

•V, I ; Ml,* 2 . (57) 

This Huang-Rhys factor together with the ratios of the 
energies Mlj, h'jj and V determines the shape of the 
absorption spectra. 

First, we consider a spectral density that is a single 
Lorentzian centered at a frequency O. Here and in the 
following, we choose Ml as the unit of energy. Conse- 
quently, we express jj and Qj in units of O and V is 
given in units of Ml. As a criterion for the agreement 
between a ZOFE spectrum and the corresponding PM 
spectrum we use the overlap of the areas of the two spec- 
tra. An overlap of 100% then means, that the two spec- 
tra are in perfect agreement. In Figure this overlap 
is plotted against the monomer-monomer interaction V 
for the case of a dimer (N = 2) for a single Lorentzian 
spectral density with Huang-Rhys factor X = 0.64 and 
width 7 = 0.25. From this plot we see, that the overlap 
has minima at roughly V = ±0.4. It increases rapidly go- 
ing from there to smaller \V\ and increases (more slowly) 
going to larger \V\. At V = 0, the case of indepen- 
dent monomers, the overlap reaches 100%, since in this 
case the ZOFE is no approximation, but gives the exact 
result (see Appendix IA 1[) . This can also be seen in Fig- 
ure[Tb, where the ZOFE (dashed line) and PM (solid line) 
monomer spectrum (V — 0) arc shown for the considered 
spectral density J(lu) displayed in the inset of the figure. 
There, the dashed and the solid line are indistinguish- 
able, since they lie exactly on top of each other, showing 
the perfect agreement between ZOFE and PM monomer 
spectrum. The monomer spectrum consists of broadened 
peaks separated by the vibrational energy Ml of the PM. 
We have chosen the zero of the energy axis to be located 
at the mean of the monomer spectrum. The values of the 
overlap in Figure QJ, converge to 100% also for large \V\ : 
showing that in the case of strong monomer- monomer 
interaction the ZOFE approximation becomes accurate 
too. In Figure [IJi the values of V, for which the overlap 
is minimal, are indicated by two dashed vertical lines. 
Further vertical lines mark the four values of V, where 
the overlap reaches 97%. The corresponding dimer ZOFE 
(dashed line) and PM (solid line) spectra for the marked 
values of V are shown in Figure[TJ>h. Figure[lJ shows the 
ZOFE and PM spectrum for V = —1.5, having an over- 
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FIG. 1: Left column: (a) Overlap of ZOFE spectrum with PM 
spectrum over V for a dimer for a single Lorentzian spectral 
density centered at Q. with X — 0.64 and width 7 = 0.25 £7. 
(b) Monomer spectrum and spectral density J(uj) in units of 
0.1/i 2 f2 (inset), (c)-(h) Corresponding dimer ZOFE spectra 
(dashed) and PM spectra (solid) for the values of V indicated 
by the dashed vertical lines in (a). Right column: same as 
left column, but with 7 = 0.5 Q. 



lap of 97% corresponding to the left most vertical line in 
Figure QJu As we can see here, this overlap value of 97% 
represents nearly perfect agreement between the spectra. 
The spectrum is much narrower than the monomer spec- 
trum (note the different scales of the absorption axes). 
This narrowing in the case of strong interaction V (usu- 
ally termed exchange narrowing) has also appeared in 
the investigation of Gaussian disorder [5(1 Ho, HH, sin- 

t vibrational modes (50j and in semi-empirical theories 
|34| . Apparently, as for the monomer, in this case of 



strong negative V the ZOFE approximation is quite accu- 
rate. Upon increasing V, one enters the intermediate in- 
teraction regime where discrepancies between ZOFE and 
PM spectra appear. At V = —0.41, where the overlap 
has a pronounced dip (see vertical line in Figure Hk), the 
agreement between the spectra is worst, shown in Fig- 
ure [TJi. However, when V is slightly changed from this 
value the agreement increases rapidly. Upon increasing 
V, at V = —0.1 the overlap reaches again 97% (see Fig- 
ure QJ) owing to the fact, that in the region of weak 
inter-monomer interaction, where the spectra are simi- 
lar to the monomer spectrum, the ZOFE approximation 
gives again accurate results. This also holds true for the 
case of positive weak interaction, as is demonstrated in 
Figure [If and as can be seen from the overlap values in 
Figure [T]a,. Increasing V further to positive intermedi- 
ate V again leads to discrepancies between the spectra 
as in the case of negative intermediate V. However, for 
positive V the largest deviation between ZOFE and PM 
spectrum at V = +0.44 (see Figure [Tg), where the over- 
lap is 88%, is not as large as the deviation at the overlap 
minimum for negative V. For strong positive interaction, 
there is again perfect agreement between ZOFE and PM 
spectra, as is shown in Fig.QjL for V = +1.5. These spec- 
tra have a strong blue shift w.r.t. the monomer spectrum 
and have become narrower again. 

For a larger 7, which leads to a faster decay of the bath 
correlation function a(r) (sec Eq. (|4"S")) ), we expect that 
the agreement between ZOFE and PM spectra becomes 
better, since for infinitely fast decay (the Markov limit) 
ZOFE is exact (see Appendix IA 2|) . That this reasoning 
is indeed correct is demonstrated in the right column of 
Fig. [TJ where 7 = 0.5, i.e. twice as large as in the left 
column. In the right column, the minimum values of the 
overlap of ZOFE and PM spectra are not as small as in 
the left column and the deviations between the spectra 
overall have become smaller compared to the left column. 

Next, we consider the dependence of the agreement 
between the spectra on the Huang-Rhys factor X (and 
thus on the coupling strength T). In Fig. f3J overlap val- 
ues and spectra are shown for a X roughly twice as large 
as in Fig. [TJ the values of all other parameters are the 
same. The stronger coupling of the PM to the electronic 
excitation can clearly be seen in the monomer spectra 
(Fig. f2p, j) by an increase of the intensity of the second 
peak (located roughly at energy 0). The overlap of ZOFE 
and PM spectra in Fig. [2^ shows a similar behaviour 
as before, but now the minimum has a very singular 
character for 7 = 0.25. The agreement between ZOFE 
and PM spectra at the overlap minimum (see Fig. 
becomes worse compared to the case of X = 0.64 in 
Fig. HJl. However, for 7 = 0.5 the largest deviation of 
the X = 1.2 spectra is smaller than the largest deviation 
of the X = 0.64 spectra. Comparing the overlap plots 
for X = 0.64 and X = 1.2 (Fig.Hk, i and Fig. Hi, i), one 
sees that the agreement between ZOFE and PM spec- 
tra in the case of X = 0.64 increases faster when going 
from the overlap minima to larger |V| than in the case 
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of X = 1.2 (for both values of 7). This is also reflected 
in the larger values of \V\ to which one has to go in the 
case of X = 1.2 compared to X = 0.64 (for the respective 
7), to achieve perfect agreement between the spectra (see 
Fig. Ht, h, k, p and Fig. [2b, h, k, p). These observations 
show that the quality of the ZOFE approximation de- 
pends on the magnitude of the inter- monomer interaction 

V relative to the magnitude of the coupling T = (TlVL) 2 X 
of the electronic excitation to the PM. 

In Figure [3l overlap values and absorption spectra for 
the same spectral densities as in Figure []] are shown, but 
now for a trimer (N = 3). Here, the plots of the over- 
lap against V show a greater number of local extrema. 
The trimer spectrum for the strong negative coupling 

V = —1.5 looks similar to the respective dimer spectrum, 
but is slightly narrower and has an additional small peak 
located roughly at an energy of +2.1, due to a second 



allowed electronic transition. In this case, ZOFE and 
PM spectrum agree perfectly. Upon increasing V, as for 
the dimer we can identify the cases of strong and weak 
negative and positive interaction V, where there is good 
agreement between the spectra and the cases of negative 
and positive intermediate V, where the ZOFE and PM 
spectra show deviations. However, for the trimer we ob- 
serve, that in the region of strong positive V to achieve 
perfect agreement between the spectra, we have to go to 
larger V, namely V ~ +2.7, compared to V « +1.5 for 
the dimer. Broadening of the spectral density by a fac- 
tor of two (see right column) leads for the trimer spectra 
basically to the same effects and to the better agreement 
between solid and dashed spectra as for the dimer. 

Our findings from the comparison of ZOFE spectra 
with exact PM spectra for the case of a spectral density 
with a single Lorentzian can be summarized as follows: 
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1) We always observe perfect agreement between ZOFE 
and PM spectra for strong and weak negative and posi- 
tive inter-monomer interaction V. 

2) There are deviations between the spectra in the inter- 
mediate V region (and clear resonance-like local minima 
in the overlap plots), that become smaller upon increas- 
ing the width 7 of the spectral density. 

3) Increasing the coupling T of the electronic excitation 
to the PM leads to a slower ascent of the agreement be- 
tween the spectra when going from intermediate |V| to 
larger |V|. 

4) For the trimer, we find more local minima in the plots 
of the overlap against V than for the dimer. 

5) To achieve perfect agreement between the spectra in 
the region of strong positive interaction V, in the case 
of the trimer V has to be larger than in the case of the 
dimer. 

We have also performed exact calculations for more 
complex spectral densities, that are the sum of several 
Lorcntzians. Here, we have observed similar trends as 
described above for the case of the single Lorentzian 
spectral density. As an example, in Fig. |4] ZOFE and 
PM monomer and dimer spectra are shown for a spec- 
tral density, that is the sum of six Lorentzians centered 
at different frequencies flj with different Huang-Rhys 
factors Xj and widths jj. The six frequencies, where 
the Lorentzians are centered, range from % = 0.23 to 
% = 1-61 and the corresponding Huang- Rhys factors 
have different values Xj — 0.4 ... 0.24. In the left col- 
umn of Fig. [4j each Lorentzian has a width jj = 0.25 ilj 
and in the right column this width is increased by a factor 
two, i.e. jj = 0.5 ilj. The calculation of a ZOFE dimer or 
trimer spectrum for the single Lorentzian spectral densi- 
ties of the Figures QUI as well as for the spectral densi- 
ties of Fig. 2] consisting of six Lorentzians, can be done 
within a few seconds on a standard PC. However, while 
the calculation of a PM dimer or trimer spectrum for 
the single Lorentzian spectral densities considered above 
also takes only a few seconds, the calculation of a PM 
dimer spectrum for the spectral densities consisting of six 
Lorentzians considered here took about six hours. This 
increase of the computational effort makes investigations 
involving hundreds of spectra like the evaluation of the 
overlap against V in Figures [T][3] very time-consuming. 
Therefore, in Fig.|4]we simply choose values for the inter- 
action V, that in the case of a single Lorentzian spectral 
density considered before are typical for the V regions 
of perfect agreement and deviations between ZOFE and 
PM spectra. We see from Fig. |4j that as for the sin- 
gle Lorentzian spectral densities also in this case of a 
more complex spectral density we obtain perfect agree- 
ment between ZOFE and PM spectra for strong and weak 
negative and positive interaction V and deviations in the 
intermediate negative and positive V regions. And also 
here, the deviations become smaller when we go from 
the left column of Fig.|4]to the right column, i.e. increase 
the width of the Lorentzians of the spectral density from 
7j = 0.25% to 7j = 0.5%. 
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FIG. 4: Left column: (a) monomer spectrum for 
spectral density (shown in inset) that is a sum of 
six Lorentzians centered at different frequencies Oj = 
0.23, 0.42, 0.57, 1.29, 1.41, 1.61 with different Huang-Rhys fac- 
tors Xj = 0.4,0.07,0.18,0.24,0.12,0.24, where the width of 
each Lorentzian is 7j = 0.25 Qj. (b)-(g) Corresponding dimer 
ZOFE spectra (dashed) and PM spectra (solid) for different 
monomer-monomer interaction V. Right column: same as 
left column, but with jj = 0.5 Qj for each Lorentzian of the 
spectral density. 



VIII. SUMMARY AND CONCLUSIONS 

In this work we presented two different methods to cal- 
culate zero temperature absorption spectra of molecular 
aggregates, suited for situations in which the electronic 
excitation of a monomer couples to a structured phonon 
environment. Both methods are based on an open sys- 
tem approach, where the full Hamiltonian is divided into 
a "system" part, an "environment" and the interaction 
between them. 

In the first method, we use the non-Markovian quan- 
tum state diffusion (NMQSD) approach, where the sys- 
tem contains only electronic degrees of freedom. To 
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tackle functional derivatives that appear in the NMQSD 
evolution equation the zeroth order functional expansion 
(ZOFE) approximation has been utilized. Although the 
NMQSD is based on a stochastic Schrodinger equation, 
we have shown that in the calculation of the zero tem- 
perature absorption spectrum the stochastic processes do 
not enter explicitly, so that only a deterministic equa- 
tion has to be solved. Since in this method the sys- 
tem part contains only electronic degrees of freedom, i.e. 
the dimension of the considered single-exciton Hamilto- 
nian is N for an N-mei, this approach can be applied 
to quite large systems with complex structured envi- 
ronments. The absorption spectrum of a 15-mer with 
a highly complex spectral density of the environmental 
modes for instance can be calculated within several min- 
utes on a standard PC. 

To calculate the absorption spectrum exactly and find 
out the accuracy of the ZOFE approximation used by the 
first method, we applied a second method, the method of 
pseudomodes (PM) described in section [VlJ Here, vibra- 
tional modes (pseudomodes) are included into the sys- 
tem part, which couple to the electronic excitation. The 
electronic excitation now is not coupled directly to an 
environment anymore, but each PM couples to its own 
Markovian bath of vibrational modes. This allows for 
a description of the problem via a Markovian quantum 
state diffusion (QSD) evolution equation, that can be 
solved numerically exact. For the practically important 
case of a spectral density of the environmental modes 
that is a sum of Lorentzians [42|, HE H3, [sU the absorp- 
tion spectrum of the first method without applying the 
ZOFE approximation is equal to the absorption spectrum 
provided by the PM method (as we show in Appendix [Cj). 
However, in the PM method the price one has to pay is 
that for each Lorcntzian in the spectral density one PM is 
included into the system part, leading to a rapid growth 
of the dimension of the system Hamiltonian as the num- 
ber N of monomers or the number of Lorentzians in the 
spectral density is increased. Thus, using the PM method 
we are at most able to calculate absorption spectra for a 
dimer with a spectral density of about six Lorentzians or 
a trimer with about five Lorentzians etc. due to limited 
computer capabilities. Here, the use of the Markovian 
QSD is advantageous over a propagation of a density ma- 
trix, having a size that is the square of the dimension of 
the used basis. In the QSD, for the calculation of the 
zero temperature absorption spectrum, only a system of 
differential equations with the size of the basis has to be 
solved. The spectrum is then found by a simple Fourier 
transformation. For our choice of the basis, the Hamilto- 
nian becomes very sparse, allowing an efficient numerical 
propagation. 

To investigate the applicability of the first method, 
that uses the ZOFE approximation, in section IVIII we 
compared the calculated spectra with spectra calculated 
using the PM method for small aggregates consisting of 
N = 2 or N = 3 monomers for spectral densities consist- 
ing of one and six Lorentzians. We always found perfect 



agreement between ZOFE and PM spectra in the cases of 
strong or weak inter-monomer interaction (that is, when 
the interaction energy is large or small compared to the 
width of the monomer spectrum). However, we observed 
deviations between the spectra in the intermediate inter- 
action regime, but these deviations become small upon 
increasing the width of the spectral density. In partic- 
ular, for some narrow regions of the interaction energy, 
where we found a resonance-like decrease of the agree- 
ment between the spectra, the deviations became rele- 
vant. But also for these values of the interaction energy, 
the agreement becomes good again when the width of 
the spectral density is increased. When going from in- 
termediate to strong inter-monomer interaction V, the 
agreement between the spectra increases. We found that 
in order to obtain the same degree of agreement for a 
stronger coupling T of the electronic excitation to the 
vibrations, V has to be chosen larger. For more com- 
plex spectral densities consisting of several Lorentzians 
we basically observed the same trends as for a single 
Lorcntzian. Let us finally stress that over the entire range 
of the inter-monomer interaction V, the NMQSD-ZOFE 
approach leads to spectra with an overlap of more than 
80-90% with the exact spectra, as can be seen in Figs. [I]- 
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The NMQSD approach with the ZOFE approxima- 
tion seems to be well suited to bridge the gap between 
the dimer case, which can be treated numerically ex- 
act for a large range of spectral densities (e.g. via the 
PM approach) and very large aggregates where semi- 
empirical approximations (like the coherent exciton scat- 
tering (CES) approximation (3li l3~il Ir33| ) lead to excellent 
agreement between theory and experiment. 

As shown in Ref. j37j], the NMQSD-ZOFE approach 
is not restricted to the calculation of absorption spectra, 
but allows also the calculation of fully time dependent 
quantities like the probability to find excitation on a cer- 
tain monomer. Then however, in contrast to the present 
situation, stochastic processes will explicitly enter into 
the calculations. 



Appendix A: Exactly solvable cases 
1. Non-interacting monomers 

Here we demonstrate that ZOFE leads to the exact 
spectrum for non-interacting monomers. To this end 
we show that the approximate zeroth order operator 
0^ n) (i,s), which appears within ZOFE (sec Eq. 
satisfies the fundamental consistency condition Eq. (|34p 
or equivalently the evolution equation Eq. (|36[) : 

We insert 0^ n \t,s) into Eq. flSBJ). For non- interacting 
monomers, i.e. for V nm = 0, we find 

dtOi n \t,a)=Q, (Al) 
since the commutators vanish and the functional deriva- 
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tive yields zero. Equation (|A1[) together with the initial 
condition Oq\s,s) = L n (Eq. (|37P) yields the constant 

0^(t,s)=L n (A2) 

without any approximation. Thus, ZOFE is exact for 
non-interacting monomers. 

2. Markovian environment 

Within the QSD approach of section IIVI we choose the 
bath correlation function a n (t — s) to be delta-like, i.e. 

a n (t-s) = 2h 2 9 n 8(t ~ s) (A3) 

which amounts to a Markovian environment. Inserting 
Eqs. and (37]) into Eq. §Q yields 

OW(t,z')=Wn. (A4) 

We see that in this Markov limit, the exact operator 
CK n '(i,z*) is independent of z* and constant w.r.t. t. 
Thus, we find that ZOFE also contains the Markov limit. 
We insert Eq. (|A4[) into Eq. (|38|) and obtain the evolution 
equation 

d t \iP(t,z*)) = ~H sys \<p(t,z*)) 

(A5) 

+ Y J {LnZ*t, n -0 n L\ l L n )\i>{t,2,*)). 

n 

which is the well-known linear Markov QSD equa- 
tion [64[. 

Appendix B: Absorption in the pseudomode 
approach 

In this section, we show how we obtain the absorp- 
tion spectrum using the approach of pseudomodes of Sec- 
tion [VI] As the initial state of the aggregate we take 

|i(t = 0)) = | 5 el>|flPM>|6) > (Bl) 

where | g e \ ) is the product of all monomer electronic 
ground states given by Eq. The state | #pm ) de- 

notes the product of the ground states of all PMs and 
| ) is the bath state in which all bath modes are in their 
ground states. Analogously to Section ITlTI the cross sec- 
tion for absorption of light can be obtained from Eq. (|53[) 
with the dipole correlation function 

M(t) = ($(t = 0)\jl-£ e - im ' h /2-£*|<l(i = 0)), (B2) 

where fl is the total dipole operator of the aggregate, 
given by Eq. (TTTJ. As in Section HH by using Eqs. p7j) 
and HH]), we obtain 

M(t)=tf ot (*(t = 0)\*(t)) (B3) 



with the normalized aggregate states 

|*(t = 0)) = |&>|6) (B4) 

and with | ipo ) given by Eq. (|55p and 

| ) = exp(-iH e t/h)\ V(t = 0) }. (B5) 

The prefactor fi tot in Eq. (|B3[) is defined in Eq. (f2"0"]) . 
Similarly as in Section liVBl we use the expansion 

!*(*)>= / ^e-^ 2 \^ t} z*))\i). (B6) 

Here, \z njp ) = exp(z„ Jp a„ Jp )| 0„j P ) denote Bargmann 
coherent bath states, and | z } = Y[ n Ylj Tip I ^«ip )■ From 
Eqs. flB4|), |B6]), and dB3|, we obtain 

M(t)=i4 ot (i> \i>(t,z* =0)). (B7) 

We may determine the state | ip(t, z* = 0) ) using 
the QSD approach analogously to Section IIVI now for 
a memory-less (Markovian) environment (see also ap- 
pendix [A2|. First, H e of Eq. (gH) is transformed to the 
interaction representation w.r.t. H env Eq. (|50[) to find 

N 

H e (t) = H sys + EE + V nj A nj {t)) (B8) 

n=l j 

with 

A nj (t) = E ~< Jp e-^ pt ~a njP (B9) 
p 

and where 

L nj =b nj . (BIO) 

From the definition of A n j(t) Eq. (|B9[) . we get the bath 
correlation function 

a nj (t-s) = (0 11^)1^.(5) 1 6) 

= Ei^pi 2e "^ jp(t_s) ' (B11) 
p 

where | ) denotes the state in which all bath modes are 
in their ground states. 

We want to solve the Schrodinger equation 

ihdtWt)) =H e (t)\V(t)). (B12) 

Inserting the expansion Eq. (|B6[) into Eq. (|B12[) we ob- 
tain 

d t \ z*) ) = ~H sys \ $(t, z*) ) + E L nj z* tinj \4>(t, r) ) 

(B13) 
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with timc-dcpcndcnt complex numbers 



Z t,nj — ^2-^1 Kn 3P Z njp e 

p 



(B14) 



From Eq. (|52j) we have a n j(t ~ s) oc S(t — s), i.e. the 
environment of the pseudomodes is Markovian. Now, we 
make use of this fact by inserting Eq. (|52|) into Eq. (|B13|) 
and obtain 



d t \Tp(t,Z*))=-~H sys \i>(t,z*)) 



Y^Lnfzl^iitX)) 



ii -J 



(B15) 



£^nAil$(*.2*)>- 



In the case z* = equation (|B15|) yields the evolution 
equation ([55)1 of Section fVI Bl i.e. 



Since the corresponding matrix is very sparse, we are able 
to calculate the absorption spectrum taking into account 
a few PMs per monomer. For the calculation of the dipolc 
correlation function Eq. we choose the initial state 
Eq. (|55p to be real-valued. On that condition and by 
using the fact that the Hamiltonian corresponding to the 
system of equations Eq. (|B19|) is symmetric (it is equal 
to its transpose) , we can calculate the value of the dipolc 
correlation function at time 2t through [65j 



M(2t) 



2 

Mtot 



(B20) 



where the star * denotes the complex conjugate. That is, 
the time we need to propagate the wavefunction numer- 
ically shortens by a factor of two. Note, that for com- 
plex initial wavefunctions the efficient calculation scheme 
Eq. (|B20[) is no longer applicable. 



fie I i>(t, r = o)) = - -H sys \ v>(t, r = o) ) 



(B16) 

In order to obtain the absorption spectrum, this equation 
is solved for the initial state \ipo)- 



a. Numerical implementation 

We express the Schrodinger equation Eq. ([56]) in a basis 
of product states 

A 

i^>=K>nin^-> ( bi? ) 

n=X j 

with vibrational eigenstates | (3 n j ) of PM j in the elec- 
tronic ground state of monomer n, i.e. the states | j3 n j ) 
satisfy 

W.njb\ l jb n j\ Pnj ) = f&lnj Pnj\ Pnj )■ (B18) 

From the Schrodinger equation Eq. (|56[) we obtain a sys- 
tem of coupled differential equations for the components 
{6% | ip(t) ) w.r.t. the states Eq. (|BT7)) 

d t (0P\$(t)) = 

~ \ ( £ ™ + E E^^ - ifvy mj )p m j ] {6%\ i>{t) ) 



m=l j 



Min...l3 n3 -l...) 



Appendix C: Equality of NMQSD spectrum and PM 
spectrum 

In this section, we show that the absorption spectrum 
obtained from the NMQSD approach is equal to the spec- 
trum obtained from the PM method for a bath correla- 
tion function consisting of a sum of exponentials as in 
Eq. P5|) . This equivalence holds true provided the pa- 
rameters (n n j,j n j, and T n j) of the pseudomode descrip- 
tion are taken from the corresponding bath correlation 
function of the NMQSD approach. 

We start with the Hamiltonian of the PM approach 
defined by Eqs. P5|) - (|5T1) . It can be written as 



N 
n=l j 



L^bnj ) + -Hvib- 



(CI) 

Here, H sys is the Hamiltonian containing only electronic 
degrees of freedom and 



N 



H vih = E E m ni b lj b nj 
n=l j 

N 

+ E E E (Kjp**ipt&3 
n=l j p 

N 

+ EE E M 

n=l j p 



KnjpO> n j p bnj 



(C2) 



"njp^njp, 



n j i J- \ On 



(Pn...p nj +1...) 



j 

N 



(B19) 



contains the pseudomodes, their Markovian environ- 
ments, and their respective couplings. We diagonalize 
iJvibj i- c - we write 



A' 



n=l j <r 



(C3) 
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with new frequencies Q) n j a and new annihilation opera- 



tors c n j a that are a linear combination 



jo — Xnjo-bnj + Cnjpa-O/njp 



of the original ones and with 

= Xnj 



(C4) 



(C5) 



where the Xnja and Cn,,^ are complex coefficients (this 
is shown in Ref. (66j for real coefficients and can be 
extended easily to the present case). With Eqs. (|C3[) 
and jnSJ) H e Eq. jUI]) has the form of H e Eqs. ([9])-(fT2]) 
and we can derive a Schrodinger equation in the reduced 
space of H sys analogously to section IIVI Transforming 
H e from Eq. (|C1[) to the interaction representation w.r.t. 
H vih yields 



A' 



H e (t) = H sys + £ {L n Bl(t) + LlB n (t)) 



n=l 



with 



and 



B n (t) = J2^/T~b nj (t) 



b nj (t)=e iH ^ n b nj e - iH ^ h . 



(C6) 



(C7) 



(C8) 



Note that Eq. (|C6[) has the same structure as Eq. ([25 
Using Eqs. (fU5|) and (fC3|) . from Eq. ([08]) we obtain 



a 

From Eqs. JC9]) and ([C7l), we get 



(C9) 



(CIO) 



Next we insert the expansion 



*(*)) 



^ P -|a| 2 | 



^(t,z*))|z> 



(Cll) 



) 



w.r.t. Bargmann coherent states | z Tl 

exp(z njrJ cl ja )\£? lja ), where |z) = U„UjU<r\ Znj*), 
into the Schrodinger equation 



ihdt\*(t)) = H e (t)\*(t)) 



(C12) 



for the total state \^(t)). This leads to an evolution 
equation 



dt\<p(t,&*)) 



H S y S 



S\ V (t,z*)} 



(C13) 



for the state | ip(t, i*) ) in the Hilbert space of the original 
"system" with Hamiltonian H sys - In Eq. (|C13p . the time- 
dependent complex numbers 



(C14) 



and the correlation functions 

p n (t -s) = (B n (t)Bi( s )) T=Q = (6\B n (t)Bl(s)\0) 



(C15) 



are used. 

The NMQSD absorption spectrum is equal to the PM 
spectrum if the corresponding time correlation functions 
M(t) of the NMQSD approach and M(t) of the PM 
approach are equal (see Eqs. (fTSj) and ([SH]) - ). From 
Eqs. ([B3|-([B5| and Eq. (jUTT]) we get 



M(t)=ri ot (Q\{gp M \(^o\ I —e- 



\<p(tX))\*), 

(C16) 



similar to Eq. ([39]) . Because of equation (|C4[) . we have 

.9PM }| 6} = | 6). (C17) 
Using Eq. (|C17|) . equation (|C16|) yields 

M(t)=/i 2 ot (^o|^,2*=0)) (C18) 
with the initial state 



y(f = 0,8*=0)) = |Vto). 



(C19) 



The time correlation function M(t) Eq. (|C18|) of 
the PM method is equal to M(t) Eq. ((40]) of the 
NMQSD approach, if the state | tp(t 7 z* = 0) } is equal 
to | ip(t, z* = 0) } for all times t. Thus, we next show 
that Eq. (|C13[) for z* = is equivalent to Eq. ([30]) for 
z* = 0. According to Eq. ([55]). we replace the func- 
tional derivative in Eq. (|30[) by the operator O^ 1 ' (t, s, z*) 
and the functional derivative in Eq. (|C13[) by an opera- 
tor 0( n \t, s,z*). Note, however, that we do not make 
the ZOFE approximation, so that our treatment remains 
exact. Following Ref. [39|, one can expand the opera- 
tors CK n )(t,s,z*) and 6^(t,s,z,*) w.r.t. the function- 
als z t * n and z t * n and obtain a hierarchy of differential 

equations for the different orders , , . . . of the 
expansion. These differential equations do not contain 
zl n (and z t * n respectively) anymore and they are the 

same for 0^(t,s,z*) and 6 {n \t,s,i*). Thus we have 



°o n> = O ( n \ 0[ n> = 6[ n> , .... Now it only remains to 
show that in Eqs. ([20]) and (]C13p a n (t - s) = (3 n (t - s). 
This we will show in the following. We consider the time 
derivative of Eq. ([C8]) 



d t b nj (t) = l -e iH ^/n [Hv . hj Kj ] e -*H vlb t/n_ (C20) 
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Inserting the definition of H vih Eq. (fC2"| into Eq. (|C20j) 
yields 

d t b nj (t) = -m nj b nj (t) K* njp a njp (t), (C21) 



with 



iH vib t/h ~ . p -iH vib t/h 



(C22) 



Taking the time derivative of Equation (|C22|) and insert- 
ing the definition of -ff v ib Eq. (|C2[) yields 

dtQ>njp\t) = ~i&njpQ'njp(t) ~ T^nj pbnjXp) • (C23) 

We integrate Eq. (fC23|) and obt am 



l -k njp [ ds e-* 5 -^ (*-)6 n3 -(a). 



(C24) 



Inserting a n j p (t) Eq. (|C24|) into Eq. (|C21j) yields 
d t b nj {t) = -id nj b nj (t) 



Z 7 
ft 



(C25) 



with A n j(t) given by Eq. (|B9[) . The A n j(t) are correlated 
with the function & n j(t — s) given by Eq. (IB1 1[) . 

Using the PM method the bath is Markovian, so that 
a n j(t — s) is proportional to a delta function Eq. (|52[). 
Inserting Eq. ^ into Eq. (|C25j) leads to 



d t b n j(t) = -ifl nj b nj (t) - "( n jb n j(t) - -A nj (t). (C26) 



Integration of Eq. (|C26[) yields 
M*) =e-^ t -^*6„ j (0) 



(C27) 



We now insert b n j (t) Eq. (|C27|) into the correlation func- 
tion 



<M*)W>t=o = (0 |M*)£j(«)l 0) (C28) 

and obtain (by separately considering the cases t < s and 
t > s) 



(bn j (t)bt j ( S )) T = =e- iQ ^ t - s ^ t - 



(C29) 



Finally, Eq. (|C29| toge ther w ith the definition of B n (t) 
from Eq. (|U7f and Eq. (|C15|) yields the result 



/3„(t - a) =rV- ffl «'N-i«IH i (C30) 



i 

which is equal to the special bath correlation function 
a n (t - s) given by Eq. flU) for the NMQSD approach. 

Thus, we have shown that the absorption spectrum 
we obtain from the PM method of section I VII is equal 
to the spectrum we get by using the NMQSD approach 
according to section IIVI 
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